module mutual_information

use tests

contains

real function mi_p(probx, proby, probxy, sizex, sizey)

	implicit none

	integer, intent(in) :: sizex, sizey
	real, intent(in), dimension(0:sizex) :: probx
	real, intent(in), dimension(0:sizey) :: proby
	real, intent(in), dimension(0:sizex, 0:sizey) :: probxy
	real :: tempreal
	integer i, j

	mi_p = 0.0

	do i = 0, sizex
		do j = 0, sizey
			tempreal = probxy(i,j)*log( probxy(i,j)/(probx(i)*proby(j)))
			if( (test_isinf(tempreal) .eqv. .false. ) .and. test_isnan(tempreal) .eqv. .false.) then
				mi_p = mi_p + tempreal
			end if
		enddo
	enddo

	mi_p = mi_p / log(2.0)

	if(mi_p < 0.0) then
		write(*,*) 'MI less than 0, check head count'
		mi_p = 0.0
	end if
	

end function mi_p



real function st3mi(probx, proby, probz, probxyz, sizex, sizey, sizez)

	implicit none

	integer, intent(in) :: sizex, sizey, sizez
	real, intent(in), dimension(0:sizex) :: probx
	real, intent(in), dimension(0:sizey) :: proby
	real, intent(in), dimension(0:sizez) :: probz
	real, intent(in), dimension(0:sizex, 0:sizey, 0:sizez) :: probxyz
	real :: tempreal
	integer i, j, k

	st3mi = 0.0

	do i = 0, sizex
		do j = 0, sizey
			do k = 0, sizez
				tempreal = probxyz(i,j,k)*log( probxyz(i,j,k)/(probx(i)*proby(j)*probz(k)) )
				if( (test_isinf(tempreal) .eqv. .false.) .and. test_isnan(tempreal) .eqv. .false.) then
					st3mi = st3mi + tempreal
				end if
			enddo
		enddo
	enddo

	st3mi = st3mi/log(2.0)

end function st3mi



real function condmi(probx, proby, probz, probxz, probyz, probxyz, sizex, sizey, sizez)

	implicit none

	integer, intent(in) :: sizex, sizey, sizez
	real, intent(in), dimension(0:sizex) :: probx
	real, intent(in), dimension(0:sizey) :: proby
	real, intent(in), dimension(0:sizez) :: probz
	real, intent(in), dimension(0:sizex, 0:sizez) :: probxz
	real, intent(in), dimension(0:sizey, 0:sizez) :: probyz
	real, intent(in), dimension(0:sizex, 0:sizey, 0:sizez) :: probxyz
	real :: tempreal
	integer i, j, k

	condmi = 0.0

	do i = 0, sizex
		do j = 0, sizey
			do k = 0, sizez
				tempreal = probxyz(i,j,k)*log( (probz(k)*probxyz(i,j,k))/(probxz(i,k)*probyz(j,k)) )
				if( (test_isinf(tempreal) .eqv. .false.) .and. test_isnan(tempreal) .eqv. .false. ) then
					condmi = condmi + tempreal
				end if
			enddo
		enddo
	enddo

	condmi = condmi/log(2.0)

end function condmi


end module
